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I. INTRODUCTION 

The prediction of three-dimensional viscous flows with large separated 
regions is an essential part -of aircraft aerodynamics. For wings with 
highly swept leading edges the flow on the suction side tends to spiral in 
the manner of a vortex parallel to the leading edge. The presence of the 
rotating flow provides lift augmentation at low supersonic speeds, up to 
the point where the flow breaks down due to viscous effects. Unfortunately, 
such viscous, vortex flows do not allow easy analysis. A classical example, 
which illustrates the nature and difficulties of these flows, is the delta 
wing problem. 

The supersonic flow around a delta wing at angle of attack with sharp 
subsonic leading edges is shown schematically in Figure 1. A conical shock 
originating from the apex envelops the wing. A free shear layer separates 
from the leading edges and rolls up into a pair of conically growing 
vortices. As the angle of attack increases, the reattachment lines of 
these main vortices on the upper surface move inboard, and secondary vor- 
tices appear under the main ones, with opposite circulation. 

Previous analytical studies to solve this flow field (see Reference 1) 
have used the leading edge suction analogy (2), linear slender wing theory 
(3), or detached flow methods (4). These studies are all fundamentally 
inviscid. Some of them assume a model with two concentrated vortices lying 
on top of the wing and make use of a Kutta condition which requires the 
flow to separate tangentially from the leading edges. Thus, the viscous 
nature of the flow is contained in these conditions. Unfortunately, all 
these methods only give approximate results. A recent approach (5) uses a 
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Figure 1. General features of the flow 
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potential flow technique along with modeling of the main vortex sheet. 
However, it does not take into account secondary separation and does not 
apply as yet to supersonic flow. Finite-difference inviscid calculations 
(6) have also been performed but they do not account for the large viscous 
effects on the leeward side of the wing. 

In the present investigation, two complementary procedures are devel- 
oped which avoid the shortcomings of the above methods by solving the com- 
plete viscous and inviscid flow field about delta wings. Moreover, solu- 
tions are obtained without the costly computing requirements of a fully 
three-dimensional, time-dependent, finite-difference technique. 

In the first approach, the flow is assumed to be conically self- 
similar. This approximation is suggested by the results of experiments for 
supersonic flows around conical bodies and wings (7,8). The resulting 
Navier-Stokes equations are solved at a given Reynolds number with a time- 
marching explicit finite-difference algorithm. A similar idea has already 
been used for cones at angles of attack (9) and is currently applied to 
delta wings with supersonic leading edges (10). These calculations capture 
the bow shock, however, and are limited by a rather restrictive geometry. 
The present method treats the shock as a sharp discontinuity and allows for 
a completely general cross-sectional shape and distribution of the finite- 
difference grid points. 

In the second approach, only the streamwise viscous derivatives are 
neglected in the steady Navier-Stokes equations. This has been called the 
"parabolic" approximation because the equations take on the parabolic 
mathematical form with respect to the streamwise direction (11). The 
solution is marched downstream from a given initial station. Previous 
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investigators have used this approach, along with an implicit, iterative 
finite-difference scheme, to compute the supersonic flow over circular 
cones at angle of attack (12,13). This paper presents a new implicit, non- 
iterative algorithm which provides better computational efficiency than the 
published techniques and is not restricted to conical shapes. 

In the following pages, a detailed description of these two procedures 
is given. Some laminar results for cones and delta wings are shown and 
comparisons are made with existing computations and available experimental 
data. 
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II. FIRST METHOD: CONICAL APPROXIMATION 


Governing Equations 

The governing equations for an unsteady three-dimensional flow without 
body forces or external heat addition can be written in nondimens ional 
strong conservation- law form in Cartesian coordinates as: 


3(E-E V ) 3(F-F V ) 9(G-G V ) 

3t + 9x + 9y + 9z ° 

where E, F, G are functions of U and E v , F v , G v are functions of 

U x > Uy, U 2 . These functions are given explicitly in Appendix A. 

Conical independent variables are introduced by the following 

transformation 


( 1 ) 

U, 


ij = /x 2 + y 2 + z 2 = xA 


b, = ^ 
1 x 


z 

c, ■ — 
1 x 


( 2 ) 


where A ■ /l + bj 2 + Cj 2 . 

The conservation- law form of Equation 1 in this coordinate system 
is (14): 


i 

3b7 jll t- b l (E - V + <F-F V )]} 
+ [ ‘ C 1 (E ' 


+ bj(F-F v ) + c 1 (G-G v )] 

(E- E v ) + (F 


E v ) + (G - G v ) ] > = 0 


(3) 



The assumption of local conical self-similarity requires that deriva- 


tives of all flow quantities be zero along rays passing through the apex 
of the wing 

9E = = _9F_ _ «v = _9G_ = = 

9 a ^ 9 a ^ 9a. 9 a ^ 9 a ^ 9a^ 

This reduces the number of independent variables to three: time and 

two space variables. The calculations can be performed on a spherical 
surface centered at the apex. The viscous effects are scaled by the Rey- 
nolds number based on the radius of this surface, which is taken as refer- 
ence length L. 

Therefore a 1 ■ 1 and Equation 3 becomes 

9 / U \ 9 r-b 1 (E-E v ) + (F-F v ) 

9t V X 3; + 9b j [ X 2 

-c 1 (E-E v )+ (G-G v )“ 

+ [(E - E v ) + .b 1 (F- F v ) + c j (G - G V )J = 0 (4) 




On the sphere aj ■ 1, it is useful to define a new set of independent 
variables by the generalized transformation 

Dj = Tij (t,b 1 ,c 1 ) 

(b 1 ,c 1 ) 


(5) 


whose Jacobian is defined as 
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The final form of the governing equations in this new coordinate system 


is 


3Uj 3F : 9G 1 

Tt~ + 3nY + 


+ H 1 = 0 


(7) 


where 


U, = 


U 


1 


9ri i 1 / 3ti i 9ti i\ 

" IT U 1 " 2^2 yi 3 b 7 + C 1 9 Cl j (E ' E v ) 

l 9n i 1 9r h 

+ S^9b7 (F -V+^^-( G -Gv) 

1 / 9? 1 9C l\ 

e -§-^ b i 9b7 +C l 


1 9 ?1 1 9 ^1 
+ — (F - F„) + ^ 7- (G - G„) 


2>jX 2 9b j 


v ' ^X 2 3cj 


H 


1 ^X 


4 [(E-Ev) + b^F-Fj) + Cl (G-G v )] 


(8a) 


(8b) 


(8c) 

(8J) 


The vectors E y , F v , G y (see Appendix A) depend on IK, U^, U z> which are 
given by 

/ a„, 3nA / 3;, 3c A 

u * * '^ b > ab7 + c i 8^]% - ab7 + c i ac7j u ^! 


(9a) 


9n, 9? 1 

u y - x 9b7 % + * ab7 % 


(9b) 


9n, 9c, 

U z - > + X — u ?i 


(9c) 


•1 '1 

The system of Equations 7 is mixed hyperbolic-parabolic in time; its 
steady state solution can be obtained with a time-dependent technique. 



8 


Grid Generation 

The domain of computation on the sphere aj = 1 is limited by the bow 
shock and the body surface. Only one half of the flow field is considered; 
the other half is completed by symmetry. 

The grid required for finite-difference calculations is shown in 
Figure 2, conically projected on the physical plane (y,z) at x = 1. 

Straight rays, making an angle a with the y axis, emanate from NJ grid 
points situated at the surface of the wing. Along each ray, within the 
distance 6 between body and shock, NK points are positioned, which are 
clustered toward the wing surface. 

The choice of the surface points and angles a is arbitrary, provided 
that they are regularly distributed. In the case of the delta wing, the 
surface points are clustered toward the wing tip. The shock standoff 
distance 6 is determined by the shock boundary condition and is time 
dependent . 

The generalized coordinates rij and are defined in such a way that 
in the computational plane (nj,^) the grid has a square shape of side unity 
with uniform spacing in both directions. Therefore, the correspondence 
between physical and computational plane is, for 1 < j < NJ and 
1 < k < NK, the following: 

Pj = (k- DArij (10a) 

^ = (j - DA?! (10b) 

and 

y(j,k) = y B (j) + s(i, j ,B)cos[ct(j ) ] (11a) 

z(j,k) = z B (j) + s(i, j ,B)sin[a(j ) ] 


(lib) 
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Figure 2. Grid distribution 
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where 


An, 


A?, = 


'1 NK - 1 ’ “’’I NJ- 1 

and s is a stretching function (15) -depending on Hj , 6 , and a stretching 
parameter 6 




(lie) 


Finally, the metrics arij/abj, ar^/acj, SCj/abj, 3^/acj are obtained from 
the relations 


anj 

atT 


3c 


= ® 


i 


1 a? i 

3 n x abj 

3c7 = ‘ ®1 


H] 

ab. 


3c, 


- 2 ) 


acj \ 

i an7 
ab. 


= 0. 


1 an, 


A abj acj acj abA 


( 12 ) 


where the derivatives abj/3n^, abj/a^j, ac^/an^* aCj/a^j are computed 
numerically with central difference operators in the regularly spaced 
computational plane (16,17). 


Numerical Solution of Equations 

Numerical Algorithm 

The governing equations are solved by a time-marching finite-difference 
technique. The computations are advanced in time, from a given initial 
condition, until a steady state is reached. 

The numerical method is the standard, unsplit, explicit MacCormack 
(1969) predictor corrector algorithm (18) which has second-order accuracy 
in both time and space. A stability condition proportional to the grid 
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spacing restricts the maximum time increment. For the present viscous cal- 
culations, this time increment is computed using the empirical formula of 
Reference 19. Still, nonlinear instabilities, due to the very severe pres- 
sure gradient at the wing tip, were found to produce oscillations in the 
numerical solution. These oscillations were suppressed by using the fourth- 
order damping scheme introduced by MacCormack and Baldwin (20). 

Each time step begins by the generation of new grid and the evaluation 
of the shock boundary condition.. The finite-difference scheme is then 
implemented at each interior grid point. Finally all other boundary condi- 
tions are calculated. 

Boundary conditions 

The flow conditions at the shock boundary are computed by a "shock- 
fitting" technique. The Rankine Hugoniot relations are used across the 
shock which is allowed to move toward its steady-state position. A similar 
method was used and described in Reference 19 for a two-dimensional shock 
in body oriented coordinates. The extension to conical shocks in general- 
ized coordinates is presented in Appendix B. Beside the flow properties, 
the shock stand-off distance and the metric coefficient Sr^/St are 
obtained. 

Along the boundary ■ 0 and ■ 1, the flow variables and the 
geometric coefficients are determined using simple reflection about the 
plane of symmetry. 

At the wall, the velocities are set to zero, the temperature is given, 
and the normal pressure gradient is assumed to be zero. This assumption is 
not justified at a sharp wing tip but the loss of accuracy is minimized by 
the fine cluster of mesh points in this region. 
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Initial conditions 

The initial shock shape is an elliptic cone whose upper generator is 
a Mach line coining from the apex and whose lower generator is determined 
from a tangent-cone approximation. The initial shock speed is zero and the 
flow conditions behind the shock are obtained from the shock jump relations. 

At the wall the temperature is known. The pressure on the leeward is 
approximated by a Prandtl Meyer expansion, and on the windward by cone 
theory . 

At interior grid points, the flow variables are determined by assuming 
linear variation between the values behind the shock and those at the wall. 
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III. SECOND METHOD: PARABOLIC APPROXIMATION 


Governing Equations 

Two independent variable transformations are again applied to the 
general Navier-Stokes equations (Equation 1). They are similar to the 
transformations used for the conical approximation but allow for nonconical 
effects. The first transformation introduces conical coordinates 

a 2 = x \ 



(13) 


The second transformation allows for a stretched grid between arbitrary 
body and shock surfaces 


and 


C 2 “ a 2 

n 2 ^ 2 ^ 2*^2 1 ^* 2 ^ ► 
? 2 = C 2 (b 2 ,c 2 ) ^ 


(14) 


9 (n 2 ,s 2 ) 

" 9(b 2 ,c 2 ) 


At this point two assumptions are made: 

1. Steady state 9/9t = 0. 

2. Viscous streamwise derivatives are negligible compared with the 
viscous normal and circumferential derivatives, that is, 9/9 C 2 = 0 in the 
viscous terms only. 
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With these assumptions, the final form of the governing equations is 


3U 2 

TT 


3E 2 9F 2 

+ 9^7 + 3^7 


3G, 


= 0 


where the unsteady term 


U. 


a 2 2 U 

~W~ 


is only retained for further reference and 


a 2 2 E 

E 2 = ST 


a 2 17 
= [V 


3n- 


3n. 


a 2 3a 2 “ b 2 3b 2 


9ri 2 

+ 3b7 (F -V + 3^ (G 


8n 2 \ 

C 2 3c 2 j (E ' 

-G v )l 


^[(- b a157- < = 2 ^) (E - E v : 


9?, 

z 

3b"! 


(F-F v ) + 


The vectors Ey, F v , G v (see Appendix A) depend on U^, U , which 
given by 


3n 2 3n 2 3n 2 

a 2 3^7 “ b 2 8b7 “ C 2 3 c7 



+ 


-b 


2 3b, 


- c 


35- 
2 3Z 


Ur 


( 15 ) 


(16a) 


(16b) 


(16c) 


(16d) 

are 


(17a) 
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1 9n 2 1 
U y = ~~ TT— U n + — r Ur 

y a 2 9b 2 n 2 a 2 3c 2 £ 2 


1 ^ n 2 1 ^2 

u z - 7^ «T" u n + — Ur 
2 a 2 Sb 2 2 a 2 3 c 2 4 


(17b) 


(17c) 


This system of equations is parabolic in the £ 2 direction. It can 
be solved as an initial value problem. 

At each station x = C 2 » the generalized coordinates n 2 and ? 2 are 
defined in such a way that the domain of computation, limited by the body 
surface, the bow shock and the plane of symmetry, is mapped into a square 
of side unity. The grid generation in the computational plane (n 2 , ? 2 ) is 
identical to the one described in the above subsection on grid generation. 
It can be noted that Equation 15 is valid for nonconical body shapes. How- 
ever, for the conical shapes considered in this paper, 3n 2 /3a 2 = 0 along 
the body surface. Therefore, the body grid points can be determined in 
terms of the coordinates b 2 and c 2 only, independent of a 2 = x. 


Importance of the Streamwise Pressure Gradient 
Previous analysis 

If the initial value problem posed in the previous section is to be 
solved by forward integration in C 2 . it is clear that no upstream 
influence can be allowed in the solution. It has been shown (21) that an 
exact representation of the streamwise pressure gradient p r of Equa- 
tion 15 causes information to be propagated upstream through the subsonic 
boundary layer close to the wall. Different remedies have been proposed 

with partial success. An obvious one is to drop altogether p r from the 

^2 

equations. Cheng (11) suggests evaluating p^ with a backward difference 
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and thus to introduce it as a source term. Some authors (12,13) have 
incorporated this idea in an iterative technique, to eventually approach 
the exact representation. Rubin and Lin (12) have also proposed a '"sub- 
layer approximation" where the term p_ for the subsonic region is 
calculated at a supersonic outside the boundary layer. 

Numerical results for each procedure are given in Reference 12 for a 
two-dimensional hypersonic leading-edge problem. Except for the approximation 

p r =0, all the methods tend to exhibit instabilities and produce what is 

'2 

known as departure solutions, with a separation-like increase in wall pres- 
sure, or an expansion-like decrease in wall pressure. Lubard and Helliwell 
(13) have performed a stability analysis of their numerical scheme when 
applied to a similar system of equations. They find that the step size 

it must be greater than some minimum to avoid departure solutions. This 

s 2 

trend was verified by their numerical experimentation. 

Present analysis 

A new way of looking at this problem is to determine the influence of 
the streamwise pressure gradient on the mathematical nature of the equations 
through an eigenvalue analysis. For this, consider the two-dimensional 
parabolized Navier-Stokes equation on a flat plate, assuming constant 
velocity. A parameter u> is introduced so that these equations are 
written as 


9E* . 9P . 9F 


9F., 


9x + 9x + 9y 3y 



(18) 
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where 



PU 


_ m 

0 


pu 2 + (Jjp 


(1 - 0))p 

II 

W 

puv 

p = 

0 


[[^T p+ f (u2 + v2) ]“J 


0 


F = 


pu 


puv . 


pv z + p 


L[t^I P + 2 (^ 2 + v 2 )]vJ 




F„ = 


V 

Re 


u, ; 


4 

3 v y 


uiiy+ 2 - vv y + 


_SL 


(y - l)M 00 2 PrJ 


(19) 


In this formulation 3P/3x is to be treated as a source term with a back- 
ward difference. The problem is to determine what proportion w of p x 
can be taken out of the source term 3P/3x and included in E* without 
causing upstream influence. The inviscid limit is considered first (Re -*• <*) 

3E* , 3F 


^ + t 5 - - 0 

3x 3y 

Except for the top x term, these are the Euler equations which can be 
written also as 

A l Q x + B lQy " ' 


(20) 


0 


( 21 ) 
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where 



V 


u 

P 

0 

0* 


V 

0 

P 

o“ 


U 


0 

pu 

0 

(D 


0 

PV 

0 

0 

Q = 

V 

V 

0 

0 

PU 

0 

B l = 

0 

. 0 

pv 

1 


A 


-a 2 u 

•» 

0 

0 

u 


-a 2 v 

0 

0 

V 


( 22 ) 


and a is the speed of sound. These equations are hyperbolic in x and 
can be integrated forward in x if the eigenvalues of (A" 1 • B 1 ) are real. 
These eigenvalues are 


X 


1,2 


X 


3,4 


V 

u 

uv ± a 


✓u 2 + m(y 2 - a~^T 
u 2 - u>a 2 


and they are real if 


id 



(23a) 

(23b) 


(24) 


where M x * u/a. 

Therefore, in the region where M x > 1, the p x term can be included 
fully in E* but it must be restricted according to Equation 24 where 
M x < 1. It is only in the incompressible limit, ^ -► 0, that the entire 
pressure gradient must be in the source term. 

Next, the viscous limit is considered and the first derivatives with 
respect to y are neglected from Equation 18. In this case 


A.Q = B„Q 

2 x x 2 yy 


( 25 ) 
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where : 


A 2 " 


0 0 


0 a) 


pu 0 


u 2 + v 2 YP . P (3u 2 + v 2 ) 


B 2 


0 0 


1 0 


0 1 


(Y - l)Prp^ 


3 <y - l)Prp. 


These equations are parabolic in the positive x direction if the eigen- 
values of (AJ 1 * B 2 ) are real and positive. The eigenvalues are given by 
the zeros of the following polynomial (assuming u / 0) 


( Re IT * ' I) ( Ee V O'! 


— X) {K/ly-uiy- 1)] 




/ _ 

n 2 co ) yM x 2 

■J M x 2 + p7p — 


One can show that they will be real and positive if 


u > 0 


“ ' 1 + (y - 1)57 ■ f * (M x ) 
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The function f*(M x ) has the property that f*(l) = 1 and f*(M x ) > 1 for 
M x > 1. So that again, if M }: > 1 the term p x can be included fully, in 
E*, but must be restricted according to Equation 28b if M x < 1. Equa- 
tion 28a forbids reverse flows. 

In the present code, to is computed at each point once the flow 
variables are known. The equation for to is 

w = of* (M x ) if of*(M x ) s lj 


to = 1 


if of*(M x ) > 1 


(29) 


where o is a safety factor. 

The source term 9P/9x has not been taken into account in this analy- 
sis. It can be evaluated using a backward difference based either on the 
local pressure gradient, or on the pressure gradient outside the subsonic 
layer at a point where * Hjj >1. However, the following section shows 
that if a backward difference based on the local pressure gradient is used, 
the source term 9P/9x can have a critical influence on the stability of 
the solution, and thus may have to be dropped. 

Linear stability analysis 

In this section, the two-dimensional parabolized Navier-Stokes equa- 
tions (Equation 18) are marched in the x direction with the Euler 
implicit scheme, and a linear stability analysis of the resulting differ- 
ence equations is performed to determine which conditions must be satisfied 
by the step size Ax to obtain relaxation solutions. These conditions on 
Ax will prohibit solutions with exponential growth caused either by numer- 
ical instabilities or departure behavior. Again, only the viscous limit of 
Equations 18 is considered, but the source term 9P/9x is now included. 
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(Recall that this term represents the explicit part of the pressure 
gradient, that is, (1 - to)p x .) This system of equations can be written as: 


, 3U 9B . F iiS 

nj 9x u 9x r V^ 3y 2 


(30) 


where 


U 


P 

pu 

pv 

2 




and E^, P^, Fy represent the Jacobians 9E*/9U, 9P/9U, 9F y /9U y . 

u y 

If the Euler implicit scheme is applied to Equation 30 and the 9P/9x 
is evaluated with a local backward difference, the difference equations are 


l£ +1 - U, 1 . u/ - U^" 1 

*“5 Ax + P U Ax = F ) 


U, 


U* +1 - 2U^ +1 + U* +1 
3+1 3 3-1 

Ay 2 


(31) 


or 



where 


4 > 


Ax 

Ay2 


and the index i refers to the x direction and the index j to the y 
direction. 

In order to obtain a relaxing solution, the eigenvalues of the asso- 
ciated amplification matrix must have a modulus less than unity. The 
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coefficient matrix (13) is obtained by replacing in Equation 32 U . 1 by 
U exp (/TTKi Ay) and lA* 1 by A^'li exp (*^Tk Ay). The eigenvalues of the 
amplification matrix are the values of ' for -which the determinant of the 
coefficient matrix vanishes, 
det |) 


e 1 1 A ? - [ Ey - 2<*>F V (cos KAy - 1)] + A (P„ - Eg) - Pyj = 

' Uy 


(33) 


Equation 33 is a polynomial of degree 8 in A. If the normal velocity 
v is neglected and u is assumed to be nonzero, it can be shown that this 
polynomial can be written as 

A 3 • (A- 1) • A,X) • #(A,X,M x ,u) = 0 (34) 


where 


8 ( A 


> X) “• ^1 + | X^A - 1 


(35) 


and 


f?(A,X,M x ,w) = Ja 2 2 _ x j + x ^ Y + — 2 2 ai( ^ - ~ 1) ] - (1- oiHy- 1) | 

* { + Pr “ Y | ~ ^ “ IH^ 2 + (1 “ 2w)A.- (1 — oo) ] 


and 


"y-2 + X 2 

Pr 

Pr + X 
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2 

Pr(y- 1)M X 2 
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(y-l)M x 2 

„ _ 4yAx 

sln^) 


puReAy 2 
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(37) 
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Equation 34 has five obvious roots 

X = 0 (triple root) 

X = 1 (38) 



If u > 0, then X > 0 and these five eigenvalues always have a modulus 
less or equal to one, thus providing unconditional (neutral) stability. 

The remaining eigenvalues are the roots of the polynomial of the third 
degree % (X, X, Mjj, <o). This equation is difficult to handle analytically. 
Therefore a numerical parametric study was performed. For discrete values 
of X, M x , co such that: 

X > 0 
M x > 0 
0 < co < 1 

A numerical procedure was used to find the real and complex roots of 
Sf (X, X, Mjj, co). This procedure is a Newton-Raphson iterative technique 
where the final iteration on each root utilizes the original polynomial 
rather than the reduced polynomial to avoid accumulated errors in the 
reduced polynomial. 

From these numerical calculations, the following conclusions can be 


drawn 






(1) 

if 

M x > 1 





hi 

< 1 for all 

X > 0 

and 0 

< co < 1 

(2) 

if 

Mjj < 1 and 

is dropped 

from Equation 30 


hi 

< 1 for all 

X > 0 

and 0 

< co < f*^) 
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(It is important to note that conclusions (1) and (2) provide an indepen- 
dent check of the results obtained in the previous subsection) 

9P 

(3) if M x < 1 and — is dropped from Equation 30 

and f*(M x ) < u < 1 then j A | < 1 only if X > X min (w,M x ) 

(4) if Mjj < 1 and — is included in Equation 30 

and 0 < u> < 1 
then | X | < 1 only if X > 

These results are consistent with those of Lubard and Helliwell (13) 
who studied the two cases o> =» 0 and oj » 1 and found that 

(1) if M x > 1 and oj = 0 or ui = 1 

| X j C 1 for all X > 0 
9P 

(2) if M x < 1 and is dropped from Equation 30 

and w = 0 

1 A j < 1 for all X > 0 
9P 

(3) if M x < 1 and is included in Equation 30 

and oj = 0 
then | X ] <1 only if 


X > 


1 ~ Mx 2 

YM X 2 


(39) 


(4) if M x < 1 and oj * 1 
then | X | <1 only if 


X 


2(1 - M x 2 ) 
^x 2 


(40) 


Figures 3—6 illustrate these conclusions. In these figures, the modu- 
lus of the largest root of the polynomial ff(X) is plotted versus the 
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Figure 3. Domain of stability for to = f*(M 
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Figure 5. Domain of stability for to = 1 
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Figure 6. Domain of stability for o> = 0 and source term 
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parameters and Log X, for fixed values of u>. If the modulus of the 

largest root is greater than one, the plotting routine sets it equal to one. 
With this procedure, the regions of instability are represented by a flat 
surface which is easy to detect. 

Figures 3 and 4 illustrate the role of the pressure gradient on stabil- 
ity. Here the parameter u> is determined by Equation 29: 

If M X <1 w - f*^) 

If Mjj > 1 u « 1 

In Figure 3 the source term 3P/3x is included and there is a region of 
instability for small X and Mjj. If the source term 3P/3x is dropped 
(Figure 4), this region of instability disappears. Figures 5 and 6 compare 
the results of the present analysis with those of Lubard and Helliwell. In 
Figure 5 the parameter w is set equal to one (completely implicit pressure 
gradient). Again, there is an unstable region at small X and M x . The 
limit of this unstable region, as determined by Lubard and Helliwell 
(Equation 40) is also shown. It is clear that both analyses agree very 
well. This is also true for the case u> * 0 (completely explicit pressure 
gradient) as can be seen in Figure 6. 

At this point, it is necessary to look at the physical meaning of the 
existence of a minimum value for the parameter X. The analysis of this 
section is a viscous analysis, therefore strictly valid only for the first 
point off the wall boundary. This point is situated at a distance Ay 
above the wall. For simplicity, the boundary condition at the wall can be 
taken as a Dirichlet boundary condition (fixed U wa ^). Then the numerical 
solution of the difference Equation 32 will generate a round-off error 
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whose dominant harmonic is likely to have a period of 4Ay. If only this 
harmonic is considered, then 



(42) 


and 


v = 2y ,Ax 
PuRe Ay 2 


(43) 


The condition 

x > \in 

is equivalent to imposing a lower bound on the marching step. In particu- 
lar, if the mesh Reynolds number. Re = (Pu Ay)/u, is taken as unity for accu- 
racy purpose, it turns out that: 


2 If > Kmin < 44 > 

This unusual stability condition has been verified experimentally by Lubard 
and Helliwell (13). Unlike Reference 13, the present analysis does not 
provide an analytical formula for the lower bound on Ax. However, this is 
not so restrictive since in a real problem, the minimum Ax has to be 
determined by trial and error. 

As a conclusion, it is important to recall the main result of the 
present analysis: the best strategy to obtain unconditional stability is to 

include only part of the pressure gradient in the normal implicit algorithm — 
namely, cop x where w is given by Equation 41 — and drop the other part 
entirely, that is, (1 - u>)p x . 
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Numerical Solution of Equations 

Numerical algorithm 

Equation 15 is solved with a finite-difference technique adapted from 
the class of completely implicit, noniterative ADI schemes introduced by 
Lindemuth and Killeen (22), Briley and MacDonald (23,24), and Beam and 
Warming (25-27). It uses the implicit approximate factorization in delta 
form of Beam and Warming (26). The choice of an implicit algorithm is justi- 
fied when the limit imposed on the marching step by the stability condition 
of an explicit method is smaller than the limit required for accuracy. This 
is the case of the delta wing at angle of attack where the gradients in the 
longitudinal marching direction are very small compared to the large normal 
gradients due to viscosity and the large lateral gradients near the tip of 
the wing. Moreover the noniterative character of the present method is 
expected to provide better efficiency than the iterative schemes of Rubin 
and Lin (12) and Lubard and Helliwell (13). 

For the governing Equation 15, written as 


8E2 3P 2 3F 2 3G 2 

3C7 + 3c7 + TnJ + 3^“ ° 


(45) 


the delta form of the algorithm for constant step size A£ 2 is 



A£ 2 /3F- SG.V- 0 . 

T+eJ \S + ^7/ + TTeJ a1 e 2 


a s p 


2 


(46) 
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where the superscript i refers to the level £ 2 = -and U 2 = a 2 2 U/S 2 

(unsteady term of Equation 15) and 

a 1 ^ - u 2 +1 - Uj 1 

and the derivatives 3 n and 3 r are approximated with central difference 

n 2 

operators. 

This algorithm has been factorized in terms of U 2 rather than E* 

because the computation of the Jacobians 3F 2 /3U 2 , 3G 2 /3U 2 is easier than 

the computation of 3F 2 /3F 2 , 3G 2 /3E*. Since the vectors E, F, G are 

homogeneous functions of degree one in U, the conservative form of the 

governing equations is maintained. 

For first-order accuracy in £ 2 » the Euler implicit scheme is used 

(9 1 = 1, 0 2 * 0). The Jacobians are evaluated at level i and A e P = A^ 1 P. 

If second-order accuracy in £ 2 . is desired, one can use the Crank-Nicolson 

scheme (0 1 = 1/2, 0 2 = 0) or the three-points backward implicit scheme 

(0j ■ 1, 0 2 = 1/2). In this case, the Jacobians should be evaluated at 

i = (1/2); this can be done through an extrapolation of levels i and i-1. 
i— 1 i— 2 

Also A g P = 2A P - A . In this study, only results obtained with the 
first-order scheme will be presented. 

The complete definition of the Jacobians is given in Appendix C. Two 
approximations are made in the computation of the viscous Jacobians. The 
coefficient of molecular viscosity is assumed to depend only on the posi- 
tion, not on the vector U. And, consistent with first-order computations, 
the cross derivative viscous terms in the (n 2 ,C 2 ) plane are neglected from 
the Jacobians. 
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In practice, algorithm 46 is implemented as follows: 



Each one-dimensional operator corresponds to a block-tridiagonal system 
of equations. In the present computations these systems are solved with a 
routine written by J. L. Steger and described in Reference 17. 

The numerical stability of the implicit portion of algorithm 46 has 
been studied by Beam and Warming for simple hyperbolic and parabolic model 
equations (27). Applied to those model equations, the Euler implicit scheme 
(6 1 = 1/2, ©2 = 0) is unconditionally stable. 

Finally, some artificial dissipation is added to the basic scheme 46. 
Fourth-order dissipation terms are added explicitly to damp eventual high- 
frequency oscillations of the solutions. These fourth-order terms are 
either identical to those used by Beam and Warming (26) and by Steger (17) 
or similar to the MacCormack damper of the conical approximation. Also, 
some second-order implicit dissipation is used. This idea, introduced to 
improve the stability of time-marching solvers (28), is used here to prevent 
departure solutions and to initiate the calculations. Its truncation error 
is consistent with a first-order Euler scheme. 
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The final algorithm is therefore (with the fourth-order explicit 
smoothing of Reference 26) 



where V and A are the conventional forward and backward difference oper- 
ators and eg and tj are the coefficients of explicit and implicit 
dissipation. 

Boundary conditions 

The conditions at the shock boundary are computed by a "shock fitting" 
technique for steady-state supersonic flows due to Thomas et al. (29). 
Details of the procedure can be found in Appendix B. The required pressure 
behind the shock is determined by an implicit one-sided integration of the 
governing Equation 15. Because this technique is not truly implicit, it 
puts a limit on the allowed integration step size A£ 2 * This limit is much 
larger than the one which would be imposed by an explicit stability condi- 
tion near the wall. However, it can be smaller than the minimum step size 
required to include the source term 9P 2 /3£ 2 ( see the linear stability 
analysis of the previous section). In this eventuality, the source term 
has to be dropped. 
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At the body surface, the change of the conservative variables 
A^U 2 j k-jjK is extrapolated from the previous value A^" and then used 

as known boundary condition in the solution of each one-dimensional normal 
operator of algorithm 46. Once the flow variables have been found at the 
interior grid points, the surface values are computed as in the conical 
approximation. The velocities are set to zero, the temperature is given, 
and the normal pressure gradient is assumed to be zero. 

The plane-of-symmetry boundary conditions are computed by reflection 
and they are imposed implicitly. 

Initial conditions 

In addition to the boundary conditions, some initial conditions are 
necessary. Ideally, the region near the apex of the wing should be com- 
puted with the full Navier-Stokes equations. In the present study the 
conical approximation described in Section II is used to generate a start- 


ing solution. The calculations are then advanced downstream to a station 


where they are compared with another conical solution and with avail- 
able experimental data. 
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IV. RESULTS AND DISCUSSION 
Test Conditions 

Laminar calculations have been performed for three test cases. A 
description of the test conditions is given in Table 1. The first case is 
a circular cone, at angle of attack, for which experimental as well as 
numerical results are available. It provides a good evaluation of each 
procedure described in the previous sections before considering the delta 
wings of cases No. 2 and No. 3. 


Table 1. Test conditions 



Test case 1 

Test case 2 

Test case 3 

Experiment 

Tracy (80) 

Monnerie and Werld (7) 

Thomann (8) 

Body shape 
Half angle or 

Circular cone 

Delta wing 

Delta wing 

sweep 

10 e 

75° 

75.7° 

Angle of attack 

24° 

10° 

9.5° 


7.95 

1.95 

3.04 

Re L 

0.42*10 6 

0.76*10 6 

10 6 

T wall /T * 

5.59 

1.13 

3 


Results from the Conical Approximation 

Test case No. 1 

For the cone calculations the mesh had 20 points along the surface and 
31 across the shock layer. The constant-? rays were chosen normal to the 
surface with a spacing of 10°. The stretching parameter 6 was set equal 
to 1.12.' The results are compared with the experimental data of Tracy (80) 
and the numerical calculations of McRae (9). 

Figure 7 shows a crosscut of the cone, the shock shapes and the tan- 
gential conical cross-flow velocity contours. Outlined is the zone of 
reverse crossflow. The agreement between experiment and computations is 
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Figure 7. Test case No. 1 — conical approximation cross-flow 

velocity contours 
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excellent, for the shock shape as well as the separation point. The surface 
pressure distributions are given in Figure 8. Also presented is the surface 
pressure computed with the conical approximation at a station situated at 
20% of the length of the cone (R.e^ = 0.84x10^). It is not the same as the 
pressure computed at Re * 0.42*10 6 . This illustrates the paradox of the 
conical approximation applied to viscous flows: the calculations remain 

Reynolds -number dependent. 

Test case No. 2 

The grid used for the delta wing calculations is presented in Figure 2. 
It has 36 points along the surface and 50 across the shock layer, with 
6 = 1.05. The numerical results are compared with the experimental data of 
Monnerie and Werle (7). Those results were obtained with fourth-order 
damping coefficients equal to 0.4 in both n and ? directions. Figure 9 
shows a crosscut of the wing, the calculated shock shape and pressure con- 
tours, along with the experimental shock position in the plane of symmetry. 
The calculated surface pressure distribution is shown in Figure 10; since 
no data are available, it is compared with Prandtl Meyer expansion for the 
leeward and with inviscid cone theory for the windward. Figure 11 shows the 
Cartesian cross-flow velocity directions immediately above the wing (the 
scale in the normal direction is twice the scale in the tangential direc- 
tion). The agreement with the experimental position of the main vortex is 
excellent. One can also see the small region of secondary separation near 
the tip. Pitot pressure measurements have also been performed in the region 
of the vortices. The pitot tube was parallel to the wing axis so that these 
measurements may be inaccurate, because of the large cross-flow velocities. 
In Figure 12, the data are compared with the computed pitot pressures based 
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Figure 8. Test case No. 1 — conical approximation surface pressure 
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Figure 9. Test case No. 2 — conical approximation — pressure contours 
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Figure 10. Test case No. 2 — surface pressure 
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on the component of velocity parallel to the wing axis. As expected, the 
comparison is only approximative. Figure 13 shows the tangential velocities 
along a row of grid points immediately above the lee side of the wing and 
gives the location of the separation and reattachment of the secondary 
vortex. The disagreement with the experimental location is believed to be 
due to a relatively coarse computational grid used. However, in view of the 
relatively large Reynolds number * the comparisons may be complicated by the 
presence of turbulence in the experiment. The conical crossflow Mach 
number contours are given in Figure 14; only a small portion of the conical 
crossflow is supersonic. Figures 15 and 16 show the streamwise conical 
velocity and temperature profiles along three constant-? rays emanating 
from the wing. Ray j • 1 is close to the windward plane of symmetry, 

Ray j « 30 goes through the main vortex, and Ray j * 36 is close to the 
leeward plane of symmetry. (A more exact definition of each of these rays 
is given in Table 2.) The inviscid portion of the flow field as well as 
the large viscous features (main vortex) are resolved properly. However, 
it should be noted that the numbers of grid points in the boundary layer 
(3 or 4) is not sufficient to give accurate shear stress and heat-transfer 
data at the wall. 


Table 2. Geometric data for 
figures 15 and 16 


Ray j= 

1 

30 

36 

y B = 

0.0685 

0 

0 

Z B = 

-.0123 

.1336 

-.01305 

a ® 

-2.65° 

150.88° 

182.65° 

Note: 

y B and Zg 

are given in the 

plane 

X - 1. 
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Figure 15. Test case No. 2 — comparison of streamwise conical 

velocity profiles 
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Figure 16. Test case No. 2 — comparison of temperature profiles 
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Test case No. 3 

In this experiment by Thoman (8), measurements were made on half a 
delta wing placed on the side wall of the wind tunnel. For the computa- 
tions, a grid similar to that of the previous case was used. The calcu- 
lated and experimental surface pressure distributions are compared in 
Figure 17. Theory and experiment agree very well in the outboard portion 
of the wing but not in the center portion. This difference is believed to 
be caused by the wind-tunnel wall boundary layer (Figure 17) which extends 
over half the wing and is interacting with the flow around the wing. 

Results for the Parabolic Approximation 

Preliminary testing of the parabolic approximation was made by calcu- 
lating the boundary- layer flow over a flat plate. The two-dimensional 
parabolized Navier-Stokes Equation 18 were solved using a simple Euler 
implicit finite-difference scheme. The conditions at the outer edge of the 
boundary layer were chosen as M,,, = 4 and Re^ ■ 1.9*10 6 where L is the 
length of the plate and the nondimensionalizing length. The wall tempera- 
ture was taken as T w * T*, and the viscosity was kept constant. The cal- 
culations were started at station x = 0.2 (assuming a trapezoidal velocity 
profile) and advanced to station x » 1. The parameter ^(Mjj) was evaluated 
by Equation 29 where the safety factor a was set equal to 0.9. The term 
3P/9x was approximated by a local backward difference. The results are 
compared with those of a standard boundary layer code (31) . Figure 18 
shows the streamwise velocity and temperature profiles. The agreement is 
very good and the slight differences in the region of higher temperature 
are believed to be due to the constant viscosity assumption. 




Figure 17. Test case No. 3 — conical approximation: pressure on 

upper surface 
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Test case No. 1 

The full three-dimensional code described in Section III was then 
applied to the cone at angle of attack of test case No. 1. The finite- 
difference grid was identical to the one used for the conical calculations. 
The solution was inarched from £ 2 * 0.2 to £ 2 * 1* Conical results at 
= 0*2 were taken as starting condition. Because the grid grows almost 
linearly with £ 2 , t * ie step size A£ 2 was chosen proportional to £ 2 - 
The ratio A£ 2 /£ 2 = 0.006 was determined experimentally by requiring 
that the "shock fitting" procedure be stable. The smoothing constants 
Eg and Ej were such that Eg * 1.04 A£ 2 and e^. = 8.33 A£ 2 . The 
parameter u) was calculated from Equation 29 with a safety factor of 
0.8. The 3P 2 /H 2 term was dropped from Equation 45. Figure 19 shows 
a crosscut of the cone and the bow shock, along with the tangential 
conical cross-flow velocity contours. The agreement with the experimental 
shock shape and separation point is again excellent. Also the velocity 
contours are almost identical to those obtained from the conical approxi- 
mation (Figure 7). The surface pressure distribution is presented in 
Figure 20, and it compares very well with experiment and calculations 
performed with the Lubard and Helliwell code. Figure 21 shows the variation 
of the shear stress with £ 2 > in planes situated 5° off the plane of sym- 
metry, on the leeward and windward of the cone. In logarithmic coordinates, 
they are compared with a straight line of slope (-1/2) , which corresponds 
to the classic boundary-layer result. The deviation of the results from a 
straight line for the leeward may be due to the presence of cross flow. 

The short oscillation at the beginning of the calculations is a transient 
phenomenon caused by the approximate nature of the starting solution. 
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Figure 20. Test case No. 1 — parabolic approximation: surface pressure 
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Figure 21. Test case No. 1 — parabolic approximation streamwise variation 

of the normal shear-stress 
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Some experimentation was done with the 3P 2 /3£ 2 term. If approxi- 
mated with a local backward difference it leads to quickly departing solu- 
tions. It was not possible to cure this problem by increasing the step 
size A S 2 since this would have made the shock fitting procedure unstable. 
With the sublayer approximation, slowly oscillating or departing solutions 
were obtained for 1 < < 2.5. For M Xg > 2.5 the results were within 

5% of those obtained with 3P 2 /3£ 2 = 0* 

Test cast No. 2 

For the delta wing, the solution was started from conical results at 
£ 2 * 0.5 and advanced to £ 2 = 1, with the same grid as iii the conical 
calculations. Again the step size was allowed to grow linearly with £ 2 * 
However, in this case a more severe restriction on A£ 2 was necessary to 
prevent instabilities in the wing tip region so that A£ 2 = 0.001 £ 2 . 

These results, apparently contradictory with the unconditional stability 
property of the implicit method, may be explained by the strong non- 
linearities in the vicinity of the tip. The smoothing coefficients were 
chosen so that e^, ■ 100 • A£ 2 (MacCormack smoothing) and e^ = 50 A£ 2 . 

The parameter to was computed from Equation 29 with a = 0.8. The 
term 9P 2 /3£ 2 was set equal to zero. The results are close to those 
obtained with the conical approximation. Figure 22 shows a crosscut of the 
wing and the bow shock, along with pressure contours. The surface pressure 
distribution is compared with the conical results in Figure 10. The curves 
are similar, differing only on the leeward by about 10-15%. Figure 23 shows 
the Cartesian c -s-flow velocity directions just above the wing (the scale 
in the normal direction is twice that in the tangential direction) . The 
position of the main vortex is predicted very well, but the region of 
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Figure 22. Test case No. 2 — parabolic approximation: pressure contours 
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Figure 23. Test case No. 2 — parabolic approximation: cross-flow 

velocity directions 
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secondary separation is somewhat smaller; this might be due to excessive 
smoothing and lack of resolution. This lack of resolution is again brought 
out in Figures 15 and 16 where the streamwise conical velocity and tempera- 
ture profiles along rays j ■ 1, j ■ 30, and j * 36 are compared with the 
conical' results. The agreement for the velocity profiles is excellent. 

The temperature profiles on the windward also agree very well. Some disa- 
greement appears on the leeward which is caused by the viscous terms not 
included in the conical approximation. The main differences however are in 
the boundary layer where the number of grid points is not sufficient for 
valid comparisons. 


Computation Times 

The results of this study were obtained on a CDC 7600 computer. The 
conical code required 3.61 * 10”** sec of computer time per step and per 
grid point. About 15 min were needed to obtain a solution for the cone and 
close to 2 hr for the delta wing. These numbers could be improved upon by 
using some of the recently developed algorithms (32,26). However, the 
standard MacCormack scheme was chosen for its reliability and ease of pro- 
gramming and because the main point was to evaluate the conical 
approximation. 

The parabolic code required 6.74 * 10“** sec of computer time per step 
and per grid point. This is to be compared with an average of 54 * 10”** sec 
for the Lubard and Helliwell code, thus providing a factor 8 improvement. 

The cone results took about 2 min of computer time and those for the delta 
wing less than 20 min. 
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V. CONCLUSIONS 

In this study, typical realistic three-dimensional flows with large 
separated regions have been calculated in a reasonable amount of computer 
time. Both conical and parabolic approximations have predicted quantita- 
tively the viscous and inviscid features of supersonic flows over cones and 
delta wings at angle of attack. Most notably determined is the location of 
the main vortex. The conical approach even produces results somewhat better 
than expected. However, the space-marching technique gives supplementary 
information about the streamwise variation of the flow variables and can be 
applied to nonconical bodies. 

Also presented in this paper was a new approach for solving the 
parabolized Navier-Stokes equations. A procedure was developed to avoid 
upstream influence and still retain streamwise pressure variations. Also, 
a new implicit noniterative finite-difference algorithm was implemented 
which provides substantial improvement in computational efficiency over 
previous techniques. The results prove the approach to be justified. 
However, a new shock fitting procedure will be required to remove the step 
limitation of the present method. It will then be possible to include the 
source term 9P 2 /9£ 2 (Equation 45) and thus retain the full pressure 
gradient Future work should also be directed toward calculating flow 

fields around nonconical bodies such as ogives and wing-body configurations. 
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VIII. APPENDIX A: GOVERNING EQUATIONS 

The fundamental equations governing the unsteady flow of a perfect gas, 
without body forces or external heat additions, can be written in conservation- 
law form for a Cartesian coordinate system as 

3(E-E V ) 3<F-F V ) 3 (G - G ) 

_ + + + = 0 

3t 3x 3y 3z 

where 
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pu * 
pv 
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pe t 
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+ + w 2 


pu 
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a xx 

puv 

*7" 

T xy 

puw 


T xz 
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a zz 

(pe t + p ) w 
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In addition, an equation of state must be specified. For a perfect 
gas, it can be written as 

p = (y - l)pe 

The Navier-Stokes expressions for the components of the shearing stress 
tensor and the heat-flux vector are 

■ & (h - 1 dlv ') 

"77 * fe (I? ‘ I dlv ? ) 

°« * fe (H ' 7 dI '' ? ) 

T * JL /IE + lv\ 

xy Re \3y 3 x/ 

T „ JL ( 2 * + lw\ 

xz Re V3z 3x/ 

P /3v 3w\ 

yz Re \3z 3y/ 

P 3T 

q x “ (y - l)M 00 2 Re L Pr 3x 

_ = P 3T 

q y (y - Dm^zr^pj- a y 

_ = P 3T 

q z (y - l)M 00 2 Re- L ?r 3z 

where the coefficient of molecular viscosity p is obtained from Sutherland's 
equation and the coefficient of thermal conductivity is computed by assuming 
a constant Prandtl number Pr * 0.72. 

These equations have been nondimensionalized as follows (the bars 


denote the dimensional quantities) 
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IX. APPENDIX B: "SHOCK-FITTING" PROCEDURES 


Conical Approximation 

The conical shock is allowed to move toward its steady-state position. 
The displacement of the shock is introduced through the time dependence of 
the shock standoff distance 6 in the plane x ■ 1. The problem is to 
express as a function of the fluid velocity at infinity and the rela- 

tive fluid velocity normal to the shock (see Figure 24). 

The local velocity of the shock is given by 


U s = -6 t n • N s 


(Bl) 


where iL denotes the inward unit normal to the shock 


N- 




ht 2 * W) + teJ 


(B2) 


| shock 

and the subscript shock refers to values along the shock in the plane x ■ 1. 
The algebraic value of the local shock velocity can be related to 6 t by 


U e 


. / 9z 


cos a - sin a) 

3C 1 / 


VR--^)' + fe) 2+ fe) 2 


(B3) 


| shock 

The vector component of the fluid velocity normal to and measured with 
respect to the moving shock is 

Vj « (V=o + U S N S ) • N g 


(B4) 


Substituting for V ro , U g , and N s , 6 fc 


can be obtained as 




Figure 24. Shock fitting notations 
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6 t - 


^) 2 + te ) 2 + fe ) 


3z 3y 

— cos a - t~ sm a 

d ?i 


(B5) 


shock 


1 

Finally, the metric coefficient 3rij/3t results from the differentiation 
of the stretching function (11c) 

Bn, 


26 


3t 


Ar 6 . 


t 2 - M}<m) 52 2 


(B6) 


where s is given by Equation 11c. From this point on, the method is 
identical to that described in Reference 19. 


Parabolic Approximation 

As the calculations proceed downstream, the position of the shock is 
computed simultaneously with the rest of the solution. The shock standoff 
distance 6 is obtained from the values at £ 2 through an Euler 
integration 

S(£ 2 + A£ 2 ) = 6(£ 2 ) + A£ 2 (B7) 


The problem is to determine the slope 6r at station The inward 

^2 2 

unit normal to the shock is given by 


+ #. 


3? 


3?. 


N s = 


v- * ft) - * m 


(B8) 
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where 


£ = 6 , 


cos a-^-sin «) + 


> 

-B H, 


) 


(B9) 


and the derivatives with respect to £ 2 are .taken along the shock. If Vj 
denotes the upstream flow velocity normal to the shock 


V, 2 = (N s • Vj 2 


(BIO) 


Substituting for V*, and N g , this equation can be solved for (the root 

such that 6r >0 is retained) 

**2 

("-fe v - - ^7 ) 

| + V l( u oo 2 - Vj 2 ) 

[fe) + te)] 


+ (— v B --^wV 

\h 2 v ~ J 
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• Vh 2 
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2 B h 2 " 2 
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dz 
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cos 


dy . 

cx - T-f— sin cx 


(Bli) 


3C 


* 2 w< *2 

The metric coefficient 3n 2 /3a 2 is obtained by differentiating Equa- 
tion 11c: 


3n. 

4 

3a. 


2B 






(B12) 


where s is given by Equation 11c. Once the new shock position is deter- 
mined, the application of a one-sided version of the finite-difference 
algorithm gives the pressure behind the shock. The rest of the flow 
variables result from the exact shock jump relations. 
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X. APPENDIX C: JACOBIANS 3F 2 /3U 2 , 3G 2 /3U 2 , and SEj/SI^ 

The -Jacob ians 3F 2 /3U 2 and 3G 2 / 3U 2 -are given by 


3F 2 Jj_f^2_ / 

3U 2 " 3U 2 I S> 2 
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Clearly, these Jacobians have an inviscid part and a viscous part: 
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2 V"' 2 /inviscid V“2/viscous 

The inviscid part can be written as a linear combination of 3E/3U, 3F/3U, 
3G/3U 
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The viscous part of the Jacobians is 
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where : 
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and (•>£ indicates derivatives with respect to these viscous 

Jacobians, the cross derivative viscous terms have been neglected and the 
coefficient of molecular viscosity has been assumed to depend only on the 
position, not on the vector U. 

Finally, the Jacobian 9E*/3U 2 is given by 
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